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Generating  and  demodulating  M-ary  CPFSK  using  the  FFT. 

Paul  H.  Moose 
Naval  Postgraduate  School 
Monterey,  CA. 

March  5,  1997 


1.  Introduction 

The  inverse  FFT  can  be  used  to  generate  M-ary,  full  response,  continuous  phase, 
frequency  shift  keying  (CPFSK)  and  the  FFT  can  be  used  to  demodulate  the  received 
signal  in  the  following  way.  We  begin  by  recognizing  that  an  inverse  DFT  creates  N  digital 
carrier  frequencies  with  amplitudes  and  phases  determined  by  the  frequency  domain  vector 
X  of  length  N.  M-ary  CPFSK  requires  that  for  each  symbol,  one  of  M=2‘’  carrier 
frequencies  (  with  initial  phase  equal  the  terminal  phase  of  the  previous  symbol)  be 
transmitted  for  each  symbol  where  q  bits  are  sent  per  symbol.  Therefore,  for  each  symbol, 
X  is  filled  with  all  zeros  except  a  complex  modulation  value  with  amplitude  one  and  phase 
determined  by  the  terminal  phase  of  the  previous  symbol  is  placed  in  the  position 
corresponding  to  the  frequency  to  be  sent  for  the  input  symbol.  M  of  the  N  possible 
frequencies  will  be  used,  but  only  one  frequency  will  be  used  on  any  given  symbol. 

The  receiver  simply  computes  the  FFT  of  the  received  symbols  and  extracts  from 
this  vector  of  length  N  a  vector  of  length  M  of  the  positions  corresponding  to  the  M 
transmit  frequencies.  The  FFT  amounts  to  implementation  of  a  quadrature  correlator  for 
each  of  the  N  frequencies.  We  select  from  those  the  M  that  could  have  possible  been  sent 
and  proceed  to  decode  M-ary  CPFSK  in  accordance  with  established  principles.  We  shall 
refer  to  this  type  of  modem  as  a  discrete  multi-tone  (DMT)  modem.  The  same  modem 
may  be  used  for  OFDM  and  M-ary  FSK  [Ref  1].  Diagrams  of  the  transmitter  and  receiver 
are  shown  in  Fig.  1 . 1  and  Fig.  1 .  2  below. 
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2.  M-arv  CPFSK 


Let  {s}  be  data  symbols  from  the  M-ary  alphabet  0,1,2, . ,M-1  and  let  the 

information  symbols  corresponding  to  these  be  given  by 


I  =  2*s-(M-1)  (2.1) 

such  that  the  {1}  are  from  the  alphabet  ±1,  ±3,  ±,5, . ±(M-1).  In  order  to  represent  a 

CPFSK  symbol,  we  begin  with  the  PAM  signal 

<l(t)  =  au(t-nT,).  (2.2) 

Here  u(t)  is  a  rectangular  pulse  of  amplitude  l/(2Tb)  and  duration  Tb,  the  symbol  interval.  * 

Let 

<f)(t;I)  =  271  TbdelJ  [  J  ^d(x)  dx  ]  (2.3) 

0 

be  the  time-varying  phase  of  the  carrier  which,  since  the  integral  of  d(t)  is  piecewise 
continuous,  is  piecewise  continuous.  The  baseband  CPFSK  transmission  waveform  is 

s(t)  =  A  exp{j  ct)(t;I) }  (2.4) 

and  consists  of  a  sequence  of  discrete  frequencies,  that  is  during  the  n*  symbol  interval, 

f„  =  [d{(})(t;I„)}/dt]/(27i)  =[  del_f  ]I„  /2  ,  n  <  t  <  (n+l)Tb.  (2.5) 


‘  A  rectangular  pulse  of  length  generates  “full-response”  CPFSK,  Other  pulse  shapes  generate  continuous  phase 
signals  referred  to  as  CPM  for  continuous  phase  modulation.  Longer  pulses  may  be  used,  rectangular  or  otherwise  and  they 
will  generate  “partial  response”  CPFSK  or  CPM. 
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The  {f  }  are  contained  in  the  alphabet  ±del_fi'2,  ±3del_f/2, . ,  ±(M-l)del_f/2.  We  see 

that  the  frequency  spacing  between  any  two  symbols  is  a  multiple  of  del_f,  the  phase  of 
the  waveform  is  piecewise  continuous  and  the  amplitude  is  constant.  The  phase  of  the 
waveform  during  the  n*  interval,  which  is  linear  for  CPFSK,  can  be  determined  by 
evaluation  of  (2.3)  to  yield 

(j)(t;I)=  e„  +  2ixh  I„  q(t-nTb)  (2  6) 


where 


n-i 

is  the  terminal  phase  at  the  end  of  the  (n-l)*  interval, 
h  =  del_P  Tb=  del_f  / 


(2.7) 


(2.8) 


is  the  modulation  index  and, 

q(t)  = 


(  0  t<0 

.  t/(2Tb)  0<t<T, 

,  '/2  Tb  <  t 


(2.9) 


is  the  integral  of  the  PAM  basic  pulse  u(t).  Thus  during  the  n*  interval  we  must  transmit 
the  frequency  ^  =  [d((i)(t;I„)}/dt]/(27r)  =  h  /(2Tb)  =  [del_f]In/2  with  initial  phase  0^ 

The  baseband  waveform  during  the  n*  interval  from  (2.4)  is 

s(t)  =  exp{j  4)(t;I) }  =  exp(j  e„ }  exp(27i:j  t }  (2. 10) 
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To  demodulate  CPFSK  we  compute  during  each  interval  the  complex  correlations 
for  each  of  the  possible  transmit  frequencies.  In  the  presence  of  AWGN  we  receive  during 
the  n'^'  interval  the  signal  with  complex  modulation  envelope 


r(t)  =  (2  Eb/  exp  {j  (j)(t;I)  }  +  n(t)  (2.11) 

where  is  the  symbol  energy  and  n(t)  is  the  complex  baseband  AWN  with  PSD  equal  N„ 
for  the  statistically  independent  real  and  imaginary  parts.  We  now  compute  the  complex 
correlations  for  each  of  the  M  possible  transmit  frequencies^.  The  output  of  the  k* 
correlator  for  the  n*  interval  is 

Zk  =  (1/Tb)  J^r(t)  exp{-  27rj  }  dt 
=  (2  Ey  exp{j  0„ }  exp{7rj  h  [I„  - 1®]  /2}  sinc{7rh  [I„  - 1<'‘>]/2  } 

+  n,  ,  k=0,l,2, . ,M-1,  (2.12) 

where  the  nj^  are  complex  gaussian  random  variables  with  statistically  independent  real  and 
imaginary  parts  that  have  zero  means  and  equal  variances  NyTb.  For  I„  =  I®  ,  the 
correlator  output  is  just 

z,=  (2Eb/Tb)'"exp{je„}  +  n,  (2.13) 

with  mean  output  (2  Tb)'"^’  exp{j  0„}  .The  mean  outputs  of  the  other  correlators 

depend  on  h  and  I„  .  However,  note  that  -I®  is  always  a  multiple  of  two,  so  that 
for  integer  values  of  h,  that  is  when  the  frequency  spacing  del_f  is  a  multiple  of  the 
symbol  rate  tb,  all  the  other  correlator  mean  outputs  are  equal  to  zero,  that  is  we  have 


^  The  real  and  imaginary  parts  of  the  complex  correlator  output  are  equal  to  the  in  phase  and  quadrature  outputs 
of  a  quadrature  correlator  at  frequency  f^  + 
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orthogonal  signaling.  We  shall  discuss  the  special  case  of  h  -  V2  ,  which  is  called  minimum 
shift  keying  (MSK)  shortly. 

2.1  Phase  states  [Ref  2] 

Consider  the  sequence  of  phases  generated  by  (2.7),  the  initial  phases  for  each 
symbol.  For  any  given  input  sequence  of  data,  they  form  a  phase  trajectory. 

If  h  is  a  rational  fraction,  h  =  m/p,  there  will  be  a  finite  number  of  phases  (  modulo  27i) 
that  can  occur.  These  are  the  phase  states  of  the  system.(  If  h  is  irrational,  there  can  be  an 
infinite  number  of  states.)  If  m  is  an  even  number,  there  will  be  p  phase  states 

(0,}  =  (  0, 7tm/p,  27tni/p,  37tm/p, . ,(p-l)7tm/p  }  (2. 14a) 

and  if  m  is  odd  there  will  be  2p  phase  states 

{0J  =  {  0,  rrm/p,  27i:ni/p,  Srcm/p,  .  .  .  .,(2p-l)nm/p  }.  (2.14b) 


2.1.1  MSK 

Consider  the  case  of  binary  CPFSK  with  h  =  V2,  the  modulation  index 
corresponding  to  MSK.  Since  h=del  ftr^ ,  we  note  that  the  bit  rate  is  twice  the  frequency 
spacing.  Here,  since  m  =  1  is  odd,  there  are  four  phase  states 

{0J  =  {  0, 7t/2,  It,  37t/2  }  ,  (2. 14c) 

the  message  alphabet  is 

{I}  =  (-1,  1}  (2.15) 

and  the  frequency  alphabet  is 
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{f}  =  { -del_f/2,  delJ72  }. 


(2.16) 


Assuming  the  initial  phase  state  is  0,  then  after  the  first  symbol,  the  phase, 

02  =  7ih  ,  can  be  7t/2  or  37t/2.  After  the  second  symbol,  the  phase  63  =  02 +'n:h  I2  can 
be  0,  or  n.  The  sequence  of  phase  states  can  be  depicted  as  a  trellis  as  shown  in  Fig.  2. 1 . 


0,  02  03  04  05  06 


Fig.  2.1  MSK  Phase  Trellis 


2.1.2  Viterbi  decoding 

The  bit  stream  for  CPFSK  with  a  finite  number  of  phase  states  can  be  optimally 
decoded  using  the  Viterbi  algorithm  [Ref  3].  The  Viterbi  algorithm  finds  the  most 
probable  path  through  the  trellis  diagram  in  an  efficient  manner  [Ref  4],  Consider  Fig.  2.1, 
the  trellis  for  MSK.  At  the  end  of  odd  numbered  symbols,  e.g.  symbol  number  three,  the 
system  will  be  in  the  next  even  numbered  state,  in  this  case  state  ©4 ,  which  can  be  either 
■k/2  or  3n/2.  If  it  is  in  state  ©4  =7r/2,  it  could  have  arrived  there  by  a  transition  from  state 
03  =  71:  with  an  I3  =  -1,  or  from  state  03  =  0  with  an  I3 =1 .  If  it  is  in  state  04=  37t/2,  it 
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could  have  arrived  there  by  a  transition  from  state  83  -  tc  with  an  I3  1 ,  or  from  state  63 
=  0  with  an  I3  =  -  1 . 

Consider  the  outputs  of  the  two  correlators  in  the  presence  of  AWGN  for  MSK, 
h=l/2.  From  (2. 12) 

zo  =  (2  Efe/  Tb)"'-  exp{j  83  }  exp{7rj  [I3-  (-1)]  M}  sinc(  -7i  [I3  -  (-l)]/4  )  +  (2.17) 

and 

z,  =  (2  Eb/  Tb)*'''  exp{j  83}  exp{7tj  [I3  -  (1)]  /4}  sinc{  -7t  [I3  -  (l)]/4  }  +  Uj,  .  (2. 18) 

Note  that  when  13=  (-1), 

E[zo]=(2Eb/Tb)"-exp{j83}  (2- 19) 

and 

E[z,  ]=-  j  (2  Eb/  Tb)''''  exp{j  83  }(2/7r)  (2.20) 

and  when  13=  (1), 

E[Zo  ]=  j  (2  Eb/  Tb)''^  exp{j  83  }(2/u)  (2.21) 


and 


E[z,  ]=  (2  Eb/  Tb)"'  exp{j  83 }.  (2.22) 

In  order  to  update  the  path  metric  for  state  84  =7t/2,  we  add  to  the  path  metric  at 
state  83  =  7t  the  Euclidean  distance  do|„  from  Zq  to  the  expected  value  of  Zq  with  83 -tc 
and  I3  =  -1 .  From  (2. 19)  the  expected  value  is  -  (2  .  We  add  to  the  path  metric 
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for  state  63  =  0  the  Euclidean  distance^  d,|o  from  z,  to  the  expected  value  of  z,  with  63  = 
0  and  I3  =  1  which  from  (2.22)  is  (2  EJ  .  We  select  the  minimum  of  the  two  total 
path  metrics  and  retain  the  path  with  the  minimum  metric,  deleting  the  other  path  and 
assign  the  minimum  metric  to  be  metric  for  state  64  =7t/2. 

Similarly,  in  order  to  update  the  path  metric  for  state  64  =3  71/2,  we  add  to  the  path 
metric  at  state  63  =  tt  the  Euclidean  distance  di|,j  from  Zj  to  the  expected  value  of 
with  03  =  71  and  I3  =  1,  which  from  (2.22)  is  -  (2  E^/  Tj''\  and  we  add  to  the  path  metric 
for  state  63  =  0  the  Euclidean  distance  do|o  from  Zq  to  the  expected  value  of  Zq  with  63  = 
0  and  I3  =  -1  which  from  (2. 19)  is  (2  E^/  Tb)  '"' ,  select  the  minimum  of  the  two  path 
metrics  and  retain  the  path  with  the  minimum  metric,  deleting  the  other  path  and  making 
the  minimum  metric,  the  metric  for  state  64  =37t/2. 

The  four  Euclidian  distances  are  illustrated  in  Fig.  2.2  for  the  last  case,  where  the 
system  actually  transitioned  from  state  63  =  0  to  state  64  =3n/2  with  an  I3  =  - 1 .  (See  Fig. 
2. 1)  In  this  case,  the  expected  value  of  Zq  from  (2. 19),  is  (2  Eb/  Tb)'"'  and  the  expected 
value  of  Zj,  from  (2.20),  is  -]{1EJ  Tb)^^(2/7t). 

The  update  procedure  following  even  numbered  symbols  intervals  follows  the 
same  pattern,  however,  the  expected  values  are  rotated  by  7t/2.  Note  that  at  each  stage  for 
MSK  there  are  only  two  phase  states  to  update,  even  though  the  system  contains  four 
phase  states.  In  fact  the  four  state  system  can  be  collapsed  to  a  two  state  system  simply  by 
multiplying  the  correlator  outputs  by  exp{jn7i:/2}.  This  reduces  the  path  memory 
requirements  for  the  soft  Viterbi  decoder  by  a  factor  of  two  but  the  computational 
complexity  remains  the  same  [Ref  3]. 


’  Selecting  the  path  with  the  minimum  Euclidean  distance  as  a  metric  is  equivalent  to  selecting  the  path  with  the 
maximum  aposteriori  probability  in  AWGN. 
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Re(z) 


3.  Frequency  spacing  and  symbol  rate 

For  CPFSK,  the  modulation  index  h  =  deLfi'fb  (see  (2.8))  is  typically,  but  not 
necessarily,  less  than  one.  Let  =  l/del_f  be  one  period  of  the  inter-carrier  spacing  del_f. 
Then 

Tb=h*T„  (3-1) 

where  is  the  symbol  interval.  Now  let 

P  =  ceil(h)  (3.2) 

so  that  PTo  s  Tb  gives  an  integer  number  of  cycles  of  del_f  with  a  length  greater  than  or 
equal  to  the  symbol  interval.  The  frequencies  (f)  needed  for  the  CPFSK  are  in  the 

alphabet  ±del_f/2,  ±3del_fr2, . ,  ±(M-l)del_fr2.  The  lowest  frequency  required  has 

period  Now  let 
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T  =  P*(  2TJ  =  N*At=  1/Af 


(3.3) 


be  the  period  of  the  lowest  frequency,  Af,  available  from  an  N  point  ifft  with  complex 
samples  clocked  out  at  sampling  frequency  f^  =  1/At.  The  frequencies  available  from  the 
ifft  are  0,  ±Af,  ±2Af,  ±3Af, ....  ±[(N/2)-l]Af  The  frequencies  required  are 

{  ±del_f/2,  ±3del_fr2, . ,  ±(M-l)del_f/2}  = 

{  ±P*Af,  ±3P*Af, . ,  ±(M-l)P*Af},  (3.4) 

that  is  P  multiples  of  the  lowest  frequency.  (Note  that  for  h<=  1,  P=l).  The  required 
symbol  length  is 


Tb  =h*T„  =  h*T/(2*P)  (3.5) 

which  requires  that  we  clock  out  . 

N,  =  h*N/(2*P)  (3.6) 

points  of  the  N  points  in  the  ifft.  For  example,  for  h=l/2,  P  =1,  del_fr2  =  Af  and  =N/4, 
that  is,  we  clock  out  1/4  of  the  points  in  the  ifft  and  the  phase  of  the  carrier  will  change  by 
7t/2,  which  is  correct  for  MSK.  A  MATLAB  program  that  generates  M-ary  CPFSK  as 
well  as  M-ary  FSK  and  OFDM  is  included  as  Appendix  I. 

3.1  Rational  modulation  index 

If  the  modulation  index  h  =  m/p  is  rational,  then 

N,  =  m*N/(2*P*p)  (3.7) 

can  always  be  made  an  integer  by  proper  choice  of  N.  If  h  is  an  integer,  then  P  =  h  and 
Nj  =  N/2.  If  h  is  less  than  one,  then  P=1  and 
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N,  =  m*N/(2*p) 


(3.8) 


which  will  be  an  integer  for  N  a  power  of  two  when  p  is  a  power  of  two  for  N>  2p.  Recall 
from  (2. 13)  and  (2. 14)  that  there  are  p  or  2p  phase  states  depending  on  whether  m  is  even 
or  odd.  Viterbi  decoders  with  states  equal  a  power  of  two  is  common  practice  for 
convolutional  code  decoders  [Ref  4],  Also  for  the  Army  s  AN/GRC-226  CPFSK  radio, 
h=l/2  with  bit  rates  of  256,512,1024,  and  2048  Kbps,  so  for  this  radio,  1/4  of  the  ifft 
points  are  used. 

The  point  of  this  discussion  is  that  for  many  practical  cases  of  CPFSK,  the  required 
number  of  discrete  time  points  in  the  CPFSK  symbol  specified  by  (3.6)  can  easily  be  made 
an  integer  for  N  a  power  of  two  as  required  by  the  radix  two  ffi  algorithm. 


3.2  Frequency  spacing  error 

In  the  event  that  (3.6)  is  not  an  integer  for  the  specified  h,  then  we  can  round  (3.6) 
to  the  nearest  integer  according  to 

Ns  =round[  h*N/(2*P)].  (3-9) 

As  illustrated  by  the  timing  diagram  in  Fig.  3 . 1,  if  we  wish  to  keep  the  symbol  rate  r^ 
constant  then  the  modulation  index  is  changed  from  h  to  h’  where 


h’  =  2*P*Ns  /N 


(3.10) 


and  the  frequency  spacing  is  changed  from  del_f  to  del_f  where 


The  error  introduced  by  this  roundoff  procedure  is  analyzed  in  [Ref  1],  It  is  inversely 
proportional  to  N  and  can  be  made  negligible  for  moderate  values  of  N  (~  32). 

Notice  that  (3.10)  guarantees  a  rational  h’=m/p  with  p  at  most  N/2,  so  there  are  at  most  N 
phase  states.  For  example,  32  phase  states  is  not  large  for  a  Viterbi  decoder. 


11^ - h 

T„^hT,  PT, 

II— H - h- 


h  P 


N^At  {n,  +  l)Ar 

11^ - 

N  ^At'  =  h'T^=T^ 


- 1 - ►  mod  index 

T  =N  At 


t 

T  ^NAt'  =  P{2T^) 


Fig.  3.1  Timing  relationships 


4.  FFT  receiver  for  CPFSK 


The  receiver  for  CPFSK  samples  the  baseband  received  signal  at  the  sample  rate  f^, 
generating  N^,  samples  for  each  received  symbol.  These  are  padded  with  zeros  to  length  N 
and  the  N  point  fft  is  computed.  The  output  values  found  in  the  frequency  cells 
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{  ±P*Af,  ±3P*Af, 


,  ±(M-l)P*Af) 


(4.1) 


are  extracted.  These  are  found  in  the  the  fft  output  vector  locations 

((P,  N-P),  (3P,  N-3P), . ( (M-l)P,  N-(M-l)P)}  (4.2) 

and  the  values  are  the  discrete  time  complex  correlator  outputs.  A  MATLAB  program 
implementing  a  fft  based  receiver  for  M-ary  CPFSK  as  well  as  M-ary  FSK  and  OFDM  is 
included  as  Appendix  II. 

4. 1  Discrete  time  complex  correlation 

Recall  from  section  2  that  in  the  presence  of  AWGN  we  receive  during  the  n* 
symbol  interval  the  signal  with  complex  modulation  envelope 


r(t)  =  (2  Tb)'"'  exp{j  0„  }  exp{2Tij  f„  t}  +  n(t)  (4.3) 


where  Et  is  the  symbol  energy  and  n(t)  is  the  complex  baseband  AWN  with  PSD  equal  N,, 
for  the  statistically  independent  real  and  imaginary  parts  The  discrete  time  output  of  the 
sampler  after  padding  with  zeros  is 


r(m)  =  (2  Eb/  Tj''  exp{j  0„ }  exp{27i:j  ^  At  m  }  +  n(m) 
0 


0^  m  ^  Nj-l 

Ns  <  m  <  N-1  (4.4) 


where 


f„  =  I„  del_f /2  =  I„  P*Af  =  I„  P/(N*At) 
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Now  consider  the  output  of  the  fft  in  the  frequency  cell  f'‘’=  f'‘^P*Af 
(for  convenience,  we  scale  by  l/N^ ), 

Z|(  =  (l/Nj)  Y,  r(m)exp(-2TT:j  m  I®P/N)  = 

(2  Eb/  exp{j  0„  }(l/NJ[{l-exp[  27rjN,  (I„-I«)P/N]}/(l-exp[  27ij  (I„-I«)P/N]}] 

+  n,  .  (4.5) 

This  should  be  compared  to  (2. 12),  the  analog  complex  correlation  result.  Note  that  in  the 
frequency  cell  corresponding  to  the  transmitted  symbol,  ,  and 

4  =  (2  Efe/  Tb)’''=  exp{j  0„ }+  n^  (4.6) 

as  in  (2.13)  for  the  analog  correlator.  However,  in  the  other  non-signaling  cells,  there  may 
be  a  slight  variation  due  to  the  discrete  nature  of  the  signal  and  the  zero  padding.  Using 
(3.6),  we  may  re-write  (4.5)  as 

z,  =  (2  Eb/  Tb)"‘  exp{j  0„  }(l/N,)[{l-exp[7tjh  (I„-I<^))]}/{l-exp[  Trjh  (I„-I®  )/NJ}] 

+  4  .  (4.5) 

First  consider  the  case  of  orthogonal  signaling  with  h  an  integer.  Since  is  a  multiple 
of  two  for  all  !„  *  I® ,  then  'E\z^  ]=  0  for  all  I„  *  I®  and  the  digital  complex  correlator 
outputs  are  orthogonal  as  are  those  for  the  analog  correlator  as  in  (2. 12).  However, 
consider  the  case  of  MSK  with  h=l/2  and  binary  CPFSK  with  I„-I®  =±  2  for  I„  . 

Then 

E[z,  ]=  (2  Eb/  Tb)"'  exp{j  0„  }(l/N,)[2/{  l-exp[-  Trj/NJ}]  ,  if  =  -1  (4.6) 

and 

E[zJ=(2Eb/Tb)"-exp{j0„}(l/N,)[2/(l-exp[7ij/NJ}l  ,  if  I„  =1  (4.7) 
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These  should  be  compared  to  (2.20)  and  (2.21)  for  the  analog  correlation  (  repeated  here 
for  convenience)* 

E[z,  ]=  -j  (2  E^/  exp{j  0„  }(2/7t)  ,  if  I„  =  -1  (2.20)* 

and 

E[zo  ]=  j  (2  Ey  Tb)'^’  exp{j  en}(2/7r)  ,  if  !„  =1  (2.21)* 

Note  that  for  0„  =  0  or  tc,  the  real  parts  of  (2.20)  and  (2.21)  are  zero  whereas  the  real 
parts  of  (4.6)  and  (4.7)  are  not.  However,  for  large  N,  we  can  approximate  (4.6)  and  (4.7) 
using  the  binomial  expansion  by 


E[z,  ]»-j  (2  Eb/  Tb)"'  exp{j  0„  }(2/7i)  , 

if  I„  =  -l 

(4.6a) 

and 

E[z„  ]«  j(2  Eb/  Tb)'-''  exp{j  0„  }(2/n)  , 

if  In=l 

(4.7a) 

which  are  identical  to  the  analog  correlator  outputs.  The  small  discrepancy  introduced  into 
the  correlator  channels  that  do  not  contain  the  frequency  of  the  current  symbol  is 
insignificant  to  the  receiver  performance  for  values  of  N  of  64  or  greater  (recall  that 
=N/4  for  MSK),  however,  it  is  important  to  use  the  exact  formulas  (4.6),  (4.7)  when 
debugging  code. 


5  Conclusions 

The  radix  two  fit  may  be  used  to  generate  and  demodulate  M-ary  CPFSK  with  a 
wide  range  of  rational  modulation  indices  h=m/p  so  long  as  p  is  a  power  of  two.  For 
ordinary  MSK,  h  =  V2.  The  N>M  point  flft  receiver  is  equivalent  to  a  bank  of  M  digital 
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complex  correlators  that  act  on  the  receive  baseband  signal.  The  M  complex  receiver 
outputs  may  be  used  as  the  soft  inputs  to  optimal  decoding  of  CPFSK  using  a  Viterbi 
decoder  or  they  may  be  demodulated  incoherently  directly  based  on  magnitude  or 
coherently  based  on  the  amplitude  of  the  real  part  as  in  a  conventional  receiver. 

The  primary  advantage  of  this  technique  is  its  flexibility.  With  simple  program 
control  changes,  the  same  DMT  modem  can  be  used  for  M-ary  CPFSK,  M-ary  FSK-PSK, 
and  OFDM  as  is  illustrated  by  the  MATLAB  code  in  the  Appendices. 
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APPENDIX  I 


%  DMT  Modulator 

%  Written  by:  Paul  H.  Moose 

%  Naval  Postgraduate  School 

%  Monterey,  CA 

% 

%  This  m-file  generates  M-ary  CPFSK,  M-ary  FSK  or  OFDM  transmit  waveforms. 
%It  is  implemented  using  the  inverse  fft  to  digitally  create  the  carrier 
%frequencies  required  for  each  symbol. 

% 

%INPUTS: 

%  type  ='mfskVcpfs'  or  'ofdm'  to  specily  the  desired  modulation 
%  ql  =  no.  of  bits  carried  by  each  of  the  Ml-ary  input  characters 

%  q2  =  no.  of  bits  carried  by  each  of  the  M2-ary  input  characters 

%  KK  =  no.  of  carriers  used  in  the  OFDM  symbols.  Set  to  0  for 
%  M-CPFSK  and  M-FSK. 

%  N  =  no.  of  points  used  in  the  iffi.  Must  be  greater  than  KK  for 
%  OFDM.  Should  be  at  least  4*M  (M=2''q)  for  M-CPFSK  or  M-FSK. 

%  rb  =  symbol  rate  for  M-FSK  and  CPFSK;  number  of  points  in  guard  interval 
%  for  OFDM. 

%  del_f  =  frequency  separation  of  carriers. 

% 

%  S  =  matrix  of  M-ary  input  characters  in  decimal  integer 
%  notation.  The  first  row  contains  the  Ml-ary  characters  that  will 

%  be  used  for  the  FSK,  CPFSK  or  OFDM.  If  there  is  a  second  row, 

%  it  contains  the  M2-ary  characters  for  phase  modulation  of  FSK-PSK. 

% 

%OUTPUTS: 

%  X  =  complex  baseband  output  sample  sequence. 

%  X  =  frequency  domain  array  of  modulation  values.  The  columns  are 
%  length  N  and  each  column  represents  a  transmission  symbol. 

%  X  is  formed  from  ifft(X). 

%  MM  =  frequency  domain  array  of  modulation  values  used  to  form  X. 

%  For  M-CPFSK  or  FSK,  the  columns  of  MM  are  of  length  M.  For  OFDM 

%  the  columns  of  MM  are  of  length  KK.  In  the  case  of  M-CPFSK  or  FSK 

%  only  one  of  the  rows  of  each  column  is  non-zero,  corresponding 

%  to  the  frequency  to  be  transmitted  for  that  symbol.  In 

%  the  case  of  OFDM,  all  of  the  rows  of  each  column  contain 

%  modulation  values  to  be  transmitted.  In  the  case  of 

%  M-FSK,  q  bits  are  transmitted  with  each  symbol,  while 

%  in  OFDM,  Q*KK  bits  are  transmitted  with  each  symbol. 
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%  MP  =  the  real  output  obtained  from  quadrature  modulation  of  x 
%  onto  a  carrier  frequency  fo  (fo  is  currently  set  to 

%  1200  in  the  program).  (MP  is  automatically  plotted  if 

%  there  are  fewer  than  20  output  symbols).  For  M-FSK  only. 

% 

%SUBROUTINES  REQD.  : 

%  freqa.m 

% 

%  USEAGE: 

%  [x,X,MM,MP]=dmtmod2(type,q  1  ,q2,KK,N,rb,del_f,S) 

%===—===================================== 


function  [x,X,lVlM,]VIP]=dmtmod2(type,q  1  ,q2,KK,N,rb,del_f,S) 


% 


%Initialize 

[aa,cc]=size(S); 


%M-FSK,  CPFSK;  Determine  number  of  cycles  of  fundamental  carrier  and  number  of 
%sample  points  to  be  used  to  account  for  the  fractional  cycle  when  the  del_f 
%is  not  an  exact  multiple  of  rb. 


if  type— 'mfsk' 

%if  rem(del_f,rb)==0  %Use  this  if  del_f  is  a  multiple  of  rb 
%  P=del_frrb; 

%  Ns=N; 

%else  %  Use  this  when  del  f  is  not  a  multiple  of  rb. 
cycles=del_f/rb, 

%P=fix(cycles)+ 1 ; 

P=ceil(cycles); 

fraction=P-cycles; 

Ns=N-round(fraction*N/P); 

Nss=N-fix(fraction*N/P); 

%end 
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%Display  the  Sampling  frequency ,  Exact  DEL_F,), 
%  Carson's  Rule  BW,  and  Bitrate. 
rbb=N*del_E(Nss*P); 
fs=Ns*rb 
DEL_F=fs*P/N 
%BW=rb+(2^q  1  - 1  )*DEL_F 
Ns 

Bitrate=ql*rb 

rbb 


%Generate  array  of  ones  properly  spaced  to  give  next  integer  number  of  cycles 
%  per  symbol  above  the  correct  number  of  cycles. 


s=S(l,:); 

s=P*s+l; 


%Now  compute  the  phase  modulation,  if  any  using  the  second  row  of  S 

if  aa==2 

M2=2^q2; 

delphi=2*pi/M2; 

phi=delphi*(S(2,;)); 

MP=exp(2*pi*j*phi); 

else 

MP=ones(l,cc); 

end 

for  n=l;cc 

]VfM(s(n),n)=MP(n); 


%  Locate  carriers  in  the  frequency  domain  array  of  N  digital 
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%  frequencies  with  carrier  number  one 

%at  frequency  -M/2,  carrier  M/2  at  frequency  -1, carrier  M/2+1 
%  at  zero  frequency  and  carrier  M  at  frequency  M/2-1. 

%  Minimum  value  for  N  is  2*M, 

X=freqa(N,MM); 


%Take  ifft  of  frequency  domain  array  producing  time  domain  array 
%of  cc  symbols  of  N  points  each  of  which  are  one  of  M  complex  sinusoids. 


%■ 


x=ifft(X); 

x=x(l;Ns,:);%Shorten  the  array  to  Ns  points  to  remove  the 
%fractional  cycle  when  del  f  is  not  multiple 
%of  rb. 


elseif  type  ==  'cpfs' 


h=del_frrb 

%P=fix(cycles)+ 1 ; 
P=ceil(h) 

%fraction=P-cycles; 

Ns=round(h*N/(2*P)); 

%Nss=N-fix(fraction*N/P), 

%end 


%Display  sampling  freq,  exact  freq  spacing,  and  exact  mod  index. 
% 

fs=Ns*rb 

DEL_F=2*fs*P/N; 

Ns; 

h=  2*P*Ns/N; 

mm=P*Ns; 

pp=N/2; 

while  rem(mm,2)=0  &  pp>l 
mm=mm/2; 

pp=pp/2; 

end 

dispC  Ns  h  mm  pp  DELF  ') 

disp([Ns  h  mm  pp  DEL_F]) 
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%Generate  array  of  ones  properly  spaced  to  give  next  integer  number  of  cycles 
%  per  symbol  above  the  correct  number  of  cycles. 


s=S(l,:); 

I=2*s-(2^ql-l) 


%Now  compute  the  cumulative  phase  modulation. 
%cumsum(I) 

cumphi=[0  rem(pi*h*cumsum(I),2*pi)] 

loc  =rem(  P*I+N,N)+1 

X=zeros(N,cc); 

forn=l;cc 

X(loc(n),  n)=exp(j  *  cumphi(n)) ; 
end 


%■ 


%Take  iffl  of  frequency  domain  array  producing  time  domain  array 
%of  cc  symbols  of  N  points  each  of  which  are  one  of  M  complex  sinusoids. 


% 


x=ifft(X); 

x=x(l:Ns,:);%Shorten  the  array  to  Ns  points  to  remove  the 
%fractional  cycle  when  del_f  is  not  multiple 
%ofrb. 


% 


elseif  t3q)e='ofdm' 
%=— =” 


fs=N*del_f 
BW=KK*del_f 
Bitrate=KK*ql  *del_f 
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%The  serial  symbol  stream  s  is  first  inverse  muxed  into  KK  streams 
%that  will  be  the  rows  of  the  matrix  S.  The  columns  of  S 
%will  be  ofdm  symbols. 

r=rem(cc,KK); 

ifr~=0 

dispC  ') 

disp(Tnput  being  truncated  by') 
disp(r) 

disp('symbols') 

end 

sof=S(l,l;cc-r); 

L=length(sof)/KK; 

Sof=reshape(sof,KK,L); 

[KK,L]=size(Sof); 


%■ 


%Modulation  values  MM  with  amplitude  one  and  one  of  2^q=M  equal 
%  phase  values  are  differentially  coded  for  each  of  the  KK 
%  carriers.  The  first  column  (symbol)  is  one  zero  phase. 

%  The  next  columns  are  differentially  coded  in  phase. 


dph=2*pi/2''ql; 


SD=cumsum(Sof)';  %  Differentially  code  the  phase  values 

MM=  exp(i*dph*SD);%  Generate  the  modulation  values. 

MM=[ones(KK,l)  MM];%  Add  the  reference  modulation  values. (Should 
%change  these  from  all  ones  in  the  future) 


% 


%  Locate  the  modulation  values  in  the  frequency  domain  array 
%  of  digital  carriers 

X=freqa(N,MM); 

0/^====_================— ======— : 
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%  Create  the  multiple  ofdm  carriers  by  executing  the  iffl  and  add  the  guard 
%  interval. 


x=iffi(X); 


ifrb=0 

x=x; 

else 

x=[x(N-rb+l:N,:);x]; 

end 

end 

%===================- 


%  Quadrature  modulate  the  baseband  symbols  onto  an  IF  carrier  frequency  fo  and 
%if  the  modulation  type  is  'mfsk'  or  'cpfsk' 

%if  type  —  'mfsk'l'cpfs' 

[rr,cc]=size(x); 


fo=800 

ko=fo*N/fs; 

ko=floor(ko)  %  The  digital  carrier  frequency  corresponding  to  fo 

MP=x(:).'; 

nn=0 ;  length(MP)- 1 ; 

MP=MP.  *exp(2*pi*i*nn*ko/N); 

MP=real(MP);  %The  real  output  signal 

%  Plot  output  automatically  for  short  inputs 
ifcc<=20 

t=0;length(nn)-l; 

t=t/fs; 

plot(t,MP) 

%end 

end 
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APPENDIX  n 


“/ofunction  [Mr,Y,theta]=dmtdmd2(type,q,KK,N,rb,del_f,y) 

% 

%  DMT  De-modulator 

%  Written  by;  Paul  H,  Moose 

%  Naval  Postgraduate  School 

%  Monterey,  CA. 

% 

%  This  m-file  demodulates  M-ary  CPFSK,  FSK  or  OFDM.  It  is  implemented 
%  using  with  an  ffl  to  which  is  equivalent  to  a  digital  correlator 
%  for  each  of  the  carriers. 

% 

%INPUTS: 

%  type  =  'mfskVcpfs'  or  'ofdm'  to  specify  the  modulation 

%  q  =  no.  of  bits  carrier  in  the  M-ary  CPFSK,  FSK  or  no.  of  bits 

%  carried  by  each  of  the  OFDM  carriers. 

%  KK  =  no.  of  carriers  used  in  each  OFDM  symbol.  Set  to  zero 
%  forM-CPFSKorFSK. 

%  N  =  no.  of  points  to  be  used  in  the  fft. 

%  rb  =  symbol  rate  for  M-CPFSK  or  FSK.  Number  of  points  in  guard  interval 

%  for  OFDM. 

%  del  f  =  frequency  separation  of  carriers 

%  y  =  matrix  of  input  time  domain  symbols.  Each  column  contains 
%  complex  baseband  samples  for  one  symbol. 

% 

%OUTPUTS; 

%  Mr  =  Matrix  of  recieved  modulation  values.  Each  column  is  a 
%  vector  containing  the  complex  modulation  values  of  the 

%  M=2^q  carriers  for  mfsk  and  cpfsk  or  the  KK  carriers  for  ofdm 

%  Y  =  Matrix  of  the  fft  of  y  after  y  has  been  zero  padded  to  N 
%  in  the  case  of  mfsk  and  cpfsk  or  after  removing  guard  intervals 

%  in  the  case  of  ofdm. 

%  theta  =  The  phase  states  for  cpfsk 

%USEAGE: 

%  function  [Mr,Y,theta]=dmtdmd2(type,q,KK,N,rb,del_f,y) 

% 

<)/„=:===„===:==============— =—==^ 


function  [Mr,Y,theta]=dmtdmd2(type,q,KK,N,rb,del_f,y) 

%============—==========—== 


%Initialize 
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[rr,cc]=size(y); 

%========== 


%M-FSK:  Extend  symbols  to  length  N  and  take  fft. 

if  type=- mfsk' 

if  rem(del_f,rb)==0 
P=del_f/rb; 

Ns=N; 

else 

cycles=del_fi'rb; 

P=fix(cycles)+1; 

fraction=P-cycles; 

Ns=N-round(fraction*N/P); 


fs=Ns*rb 

DEL_F=fs*P/N 

ye=[y;zeros(N-rr,cc)] ; 
Y=ffi(ye); 


%  Extract  KL  digital  carriers  from  the  frequency  domain  array  of  N  digital 
%  carrier  frequencies  and  place  in  the  received  array  R. 

%  Of  the  KL  freqs  in  R,  2^q  are  the  matched  filter  outputs  of  the  mfsk  signal. 
P 

KL=P*(2^q-l)+l; 

K=floor(KL/2) 

R=ifreqa(K,Y); 


%  Sample  the  filter  outputs  to  obtain  the  2^q  modulation  values  for  each 
%symbol. 

[rr,cc]=size(R); 

Mr=R(l:P:2^q*P,;); 
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elseif  type  =='cpfs' 


h=del_f/rb; 

P=ceil(h); 

Ns=round(h*N/(2*P)); 

fs=Ns*rb; 

DEL_F=2*fs*P/N; 

h=2*P*Ns/N; 

mm=P*Ns; 

pp=N/2; 

while  rem(mm,2)==0  &  pp>l 
mm=mni/2; 

pp=pp/2; 

end 

dispC  Ns  h  mm  pp  DEL  F  ') 

disp([Ns  h  mm  pp  DEL_F]) 

%  Now  compute  the  phase  states 

if  rem(mm,2)  =  0 


theta  =  pi*h.*(0:pp-l); 
else 

theta  =  pi*h.*(0;(2*pp-l)); 
end 


Y=ffi(y,N); 

Mr=[Y(N-(2^q-l)*P+l;2*P;N-P+l,:);  Y(P+l:2*P:(2^q-l)*P+l,;)]; 

%===========—==“===========—====—==== 

elseif  type=='ofdm' 

%The  precursor  is  removed  from  the  input  symbols  and  the  fft  is 
%computed 

y=y(rb+l  :rb+N,:); 

Y=ffl(y); 
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%  The  modulation  values  of  the  KK  digital  carriers  are  extracted  and 
%  placed  in  the  columns  of  array  R. 

K=floor(KK/2); 

R=iffeqa(K,Y); 

R=R(1;KK,:); 

%============================================== 


%Differentially  decode  the  symbols  in  the  time  domain. 
%  [rr,cc]=size(R); 

for  l=l:cc-l 

Mr(;,l)=R(:,l+l).*conj(R(:,l)); 
end 


end 
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